﻿using System;
using System.Collections.Generic;
using Allegro.MathInterface;
//using Cosmologica.MathInterface;

//using System.Linq;
//using System.Text;

namespace Allegro.Integrator
{
    public abstract class BaseIntegrator
    {
        double _s;
        private double[] _c = null;
        private double[] _d = null;

        protected double EPS = 1.0e-6;

        /// <summary>
        /// Tolerance
        /// </summary>
        public double Tolerance
        {
            get{ return EPS;}
            set{ EPS = value;}
        }

        /// <summary>
        /// This routine computes the nth stage of refinement of an extended trapezoidal rule. f is input as a delegate
        /// to the function to be integrated between limits a and b, also input. When called with n = 1, the routine returns
        /// the crudest estimate of Integral[f[x],{x,a,b}]. Subsequent calls with n=2,3,... (in that sequential order) will improve
        /// the accuracy by adding 2^n-2 additional interior points.
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="n"></param>
        /// <returns></returns>
        protected double Trapezoid(Func<double, double> f, double a, double b, int n)
        {
            double x, tnm, sum, del;
            int it, j;
            
            if(n == 1)
            {
                _s = 0.5*(b - a)*(f(a) + f(b));
            }
            else
            {
                for(it=1, j=1; j<n-1;j++)
                {
                    it <<= 1;
                }
                tnm = it;
                del = (b - a)/tnm; // This is the spacing of the points to be added
                x = a + 0.5*del;
                for (sum = 0.0, j = 1; j <= it; j++, x += del)
                    sum += f(x);
                _s = 0.5*(_s + (b - a) * sum / tnm); // This replaces _s by its refined value
            }
            return _s;
        }

        /// <summary>
        /// Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y, and
        /// an error estimate dy. If P(x) is the polynomial of degree N-1, such that P(Xa_i)= ya_i, i = 1,..,n, then 
        /// the returned value y = P(x)
        /// </summary>
        /// <param name="xa"></param>
        /// <param name="ya"></param>
        /// <param name="kStart">Start index of xa and ya</param>
        /// <param name="n"></param>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="dy"></param>
        protected void PolynomialInterpolation(double[] xa, double[] ya, int kStart, int n, double x, ref double y, ref double dy)
        {
            int m, ns = 1;
            double den, dif, dift, ho, hp, w;
            if(_c == null)
                _c = new double[n];
            if(_d == null)
                _d = new double[n];

            dif = Math.Abs(x - xa[kStart]);

            for(int i = 0; i < n; i++) // Here we find the index ns of the closest table entry
            {
                int j = kStart + i;
                dift = Math.Abs(x - xa[j]);
                if(dift < dif)
                {
                    ns = j;
                    dif = dift;
                }
                _c[j] = ya[j]; // and initialize the tableau of c's and d's
                _d[j] = ya[j];
            }
            y = ya[ns]; // Initial approx. to y
            ns--;
            for(m = 1; m < n; m++)
            {
                // for each column of the tableau, we loop over the current c's and d's and update them
                for(int i = 0; i < n - m; i++)
                {
                    int j = kStart + i;
                    ho = xa[j] - x;
                    hp = xa[j + m] - x;
                    w = _c[j + 1] - _d[j];
                    den = ho - hp;
                    if(den == 0.0)
                        throw new ApplicationException("Error in TrapezoidalIntegrator.PolynomialIntegration");

                    den = w/den;
                    _d[j] = hp*den;
                    _c[j] = ho*den;
                }
                if(2.0 * ns < (n-m))
                {
                    dy = _c[ns + 1];
                }
                else
                {
                    dy = _d[ns];
                    ns--;
                }
                y += dy;
                // After each coumn in the tableau is completed, we decide which correctin, _c or _d, we want
                // to add to our accumulating value of y, i.e. which path to take through the tableau - forking up or down. 
                // We  do this in such a way as to take the most "straight line" route through the tableau to its apex,
                // updating ns accordingly to keep track of where we are. This route keeps the partial approximations centered (insofar
                // as possible) on the target x. The last dy added is thus the error indication.
            }
        }
    }

    public class TrapezoidalIntegrator : BaseIntegrator, IIterativeFunctionIntegrator
    {
        const double EPS = 1.0e-5;
        const int JMAX = 20;

        public TrapezoidalIntegrator()
        {
            
        }

        public double Integrate(Func<double, double> f, double a, double b)
        {
            double s, olds = 0.0;
            for (int j = 1; j <= JMAX; j++)
            {
                s = Trapezoid(f, a, b, j);
                if(j > 5) // Avoid spurious early convergence
                {
                    if(Math.Abs(s-olds) < EPS * Math.Abs(olds) || (s == 0.0 && olds == 0.0))
                        return s;
                }
                olds = s;
            }
            throw new ApplicationException("Too many stpes in routine TrapezoidalIntegrator.Integrate");
        }
    }
}
